侠之大者，为老数据接盘
是时候把生信技能树的舞台交给后辈！昨天我分享了：生信初学者不要拿一套垃圾测序数据入门，号召大家也都分享一下自己踩过的坑，马上收到了一个超级有意思的投稿！
写在前面
在几乎所有模式植物转录组测序技术都做烂的今天，始终有一些植物因为种种原因鲜有问津。例如小麦，就像是开在奢侈品商场的黄焖鸡：有钱的人未必瞧得上，没钱的也压根就不会去光顾。不过这么多年过去了，总有那么几个祖传数据躺在NCBI的数据库里供人挖掘，比如今天这个PRJNA293629。
当RNA-seq遇上初学者
严格来讲这其实已经是我第三次挖掘这个数据。
说到这可能大家就纳闷了，这数据挖掘又不是美国大选的邮寄选票，您跑三遍还能当美国总统还是怎么着。
诚然，我这次挖的数据肯定是没有2020年魔幻，但这个数据结构也是相当的奇葩，就像是范伟遇到了赵本山，硬生生的给忽悠瘸了……
1米6还是1米7？
我们暂且跳过从sra下载到sra2fastq这一无趣的过程，相信各位都是老"省心"人。可以参考：使用ebi数据库直接下载fastq测序数据 , 需要自行配置好软件。这类工作交给服务器自己去跑就好，希望服务器不要不识抬举。
多年以后，面对shell，我仍然会回想起在NCBI严肃地敲下accession number的那个遥远的下午。为什么呢？因为这数据太有精神了，NCBI的统计结果就像是范伟的身高，read1 是125 bp，而read2 是75 bp。好家伙，这年头连NCBI的服务器都能是用花呗分期购买了吧，这数据显示不全是因为花呗没还完么。
头一次见双端测序数据还能不等长的，难道是本山大叔也开始搞生信了？赶紧翻原文method部分，清清楚楚地写着125 bp的双端测序。没想到NCBI也开始搞fake news。
按捺住rm -rf删库跑路的心情下载了fastqc的结果。从质控的结果来看，Read_1序列中0 bp和100 bp位置起碱基比例发生波动而Read_2的所有75 bp序列中碱基比例一直稳定，表明双端序列可能被错误地切割了。换句话说，100 bp双端测序被硬生生地切成了125 bp + 75 bp。
我觉得作者还是不够自信，你这应该叫200-bp双端测序，read1 200 bp，read2 0 bp。
既然知道问题出在错误地切割了PE序列，作为资深打工人，我相信你们可以很轻松地用ctrl C+V整理好这一亿多条序列。但既然我们有服务器，我这里简单说一下在linux下实现的方法。首先你需要安装seqkit和一杯手冲咖啡。我没办法教你怎么泡咖啡，但安装seqkit还是很容易的，你只需要输入：
#安装seqkit
conda install -c bioconda seqkit
之后不要加--split-3你就能直接获得200 bp双端测序数据（误），之后再使用seqkit的subseq命令分别切割数据的1-100和101-200 bp保存到read1和read2文件中 。最后不要忘记用gzip给文件减个肥，否则小心管理员给你五连鞭。
#sra2fastq，不加入--split-3，这样整个序列不会被自动分割
ls *.sra |xargs -I [] echo 'fastq-dump --gzip [] &' > fastq-dump.sh
cat fastq-dump.sh
bash fastq-dump.sh
#按照1:100和101：200分割数据,之后
for i in *.fastq.gz;
do
zcat ${i}|seqkit subseq -r 1:100 > ${i%%.fastq.gz}_R1.fq;
zcat ${i}|seqkit subseq -r 101:200 > ${i%%.fastq.gz}_R2.fq;
gzip ${i%%.fastq.gz}_R1.fq ${i%%.fastq.gz}_R2.fq;
done
作为一个具有良好自我管理意识的批处理命令，最后还应该有一个rm ${i}命令来删掉临时文件的。但如果之前的命令没有正确运行，而后面的rm又恰好完美地完成了工作，这种事情大家都没遇到过对吧，对吧。反正我现在再也不敢乱加了。
上图就是修正后的质控结果（图（1））和修正前后的hisat2比对成功率（图（2））了。图（2）中修正后样本数据的比对成功率明显高于修正前的（50%→95%），进一步证实了NCBI上的数据是有问题的。从图（2）可以看出双端等长了，双脚离地了，hisat2比对率重新占领高地了。另外，图（2）修正后的结果里SRR2306552和SRR2306530这两个文件比对率还是偏低（70%左右），通过查询NCBI可以知道这两个样本是小麦在盐胁迫下48 h后的RNA-seq结果。怎么说呢，可能实验设计上还需要考虑清楚，48小时后，两个品种明显都要死了……如果要分析，最好还是去掉最后两个时间点吧（这两个样本的基因差异表达结果大概率会被富集到细胞程序死亡相关的通路上去)。这时候就体现出Jimmy老师说的随时质控的重要性了。mutiqc一下不会吃亏也不会上当，说不定还能帮你避开大坑。年轻人，还是得 讲武德 随时质控啊。
写在最后
接下来就是计数流程了。试了几种方法，调了几种参数，最后算是"解决"了坑爹的multimapping问题了（至少summary上来看都比对上去了）。果然，还是要多看看各种工具官方提供的manual。异源六倍体伤不起呀，说到这，我又想去泡一杯咖啡了。
2018年，IGWSC发布了中国春小麦参考基因组序列，2020年，多个小麦基因组重测序项目也即将或已经公布。小麦的春天或许要来了，幸运的是也总有人正年轻着。奔涌吧， 打工人 后浪！
最好的入门方式
如果你也想开启自己的生物信息学数据处理生涯，但是自学起来困难重重，还等什么呢，赶快行动起来吧！参加我们生信技能树官方举办的学习班：
数据挖掘学习班第8期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。 生信爆款入门-第10期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。
生信技能树的粉丝都知道我们有一个全国巡讲的良心学习班，口碑爆棚，生物信息学入门省心省时省力！先看看大家的反馈吧：